Heterogeneous Diffusion in Highly Supercooled Liquids 
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The diffusivity of tagged particles is demonstrated to be very heterogeneous on time scales 
comparable to or shorter than the a relaxation time r a (= the stress relaxation time) in a highly su- 
percooled liquid via 3D molecular dynamics simulation. The particle motions in the relatively active 
regions dominantly contribute to the mean square displacement, giving rise to a diffusion constant 
systematically larger than the Einstein- Stokes value. The van Hove self-correlation function G 3 (r,t) 
is shown to have a long distance tail which can be scaled in terms of r/t 1 ^ 2 for t < 3r a . Its presence 
indicates heterogeneous diffusion in the active regions. However, the diffusion process eventually 
becomes homogeneous on time scales longer than the life time of the heterogeneity structure (~ 3r a ). 

PACS numbers: 64.70Pf, 66.10.Cb, 61.43Fs 



In a wide range of liquid states the Einstein- Stokes 
relation Dna/ksT = const, has been successfully ap- 
plied between the translational diffusion constant D of a 
tagged particle and the viscosity rj even when the tagged 
particle diameter a is of the same order as that of sol- 
vent molecules. However, this relation is systematically 
violated in fragile supercooled liquids JlJ-g]. The diffu- 
sion process in supercooled liquids thus remains not well 
understood. In particular, Sillescu et al. observed the 
power law behavior D oc r\~ v with v = 0.75 at low tem- 
peratures pi. Furthermore, Ediger et al. found that 
smaller probe particles exhibit a more pronounced in- 
crease of Drj/T tx D/Dse with lowering T 0, where 
Dse ~ kBT/2irrja is the Einstein-Stokes diffusion con- 
stant. In such experiments the viscosity changes over 
12 decades with lowering T, while the ratio D/Dse in- 
creases from of order 1 up to order 10 2 ~ 10 3 . In molec- 
ular dynamics simulations, on the other hand, the same 
tendency has been detected in a three dimensional (3D) 
binary mixture with N = 500 particles and in a two 
dimensional (2D) binary mixture with N = 1024 S. In 
our recent 3D simulation with N = 10 4 Q, r\ and D have 
both varied over 4 decades and the power law behavior 
D oc 7y~ a75 has been observed. Many authors have at- 
tributed the origin of the breakdown to heterogeneous co- 
existence of relatively active and inactive regions, among 
which the local diffusion constant is expected to vary sig- 
nificantly The aim of this paper is to numeri- 
cally demonstrate that the diffusivity of the particles is 
indeed very heterogeneous on time scales shorter than 
the structural or a relaxation time but becomes homo- 
geneous on time scales much longer than r a . 

A number of recent MD simulations have detected dy- 
namic heterogeneities in supercooled model binary mix- 
tures 10-13 to confirm a picture of cooperatively rear- 
ranging regions p4|] . That is, rearrangements of particle 
configurations in glassy materials are cooperative, involv- 
ing many molecules, owing to configuration restrictions. 
In particular, we have examined bond breakage processes 



among adjacent particle pairs and found that the broken 
bonds in an appropriate time interval (~ r Q ) are very 
analogous to the critical fluctuations in Ising spin systems 
with their structure factor being excellently fitted to the 
Ornstein-Zernike form The correlation length £ 

thus obtained increases up to the system size and satis- 
fies the dynamic scaling law, r a ~ £ z , with z = 4 in 2D 
and z — 2 in 3D. The heterogeneity structure in the bond 
breakage is essentially the same as that in jump motions 
of particles from cages or that in the local diffusivity, as 
will be discussed below. 

Much attention has recently been paid to the mode 
coupling theory |l^] . It is a self-consistent scheme for 
the density time correlation function and describes onset 
of glassy slowing down or slow structural relaxations con- 
siderably above T g . However, the mode coupling theory 
predicts no long range correlations. 

Our 3D binary mixture is composed of two atomic 
species, 1 and 2, with N\ — N 2 = 5000 particles 
with the system linear dimension L = V 1 / 3 being fixed 
at 23.2<7i. They interact via the soft-core potentials 
v a b(r) = e(a ab /r) 12 with o ah = (a a + cr 6 )/2, where r is 
the distance between two particles and a, 6 = 1,2 p6[ . 
The interaction is truncated at r = ia\. The mass 
ratio is m^/rni = 2. The size ratio is o^/ci = 1.2, 
which prevents crystallization at least in our computa- 
tion times. We fix the particle density at a very high 
value of (Ni + N 2 )/V = 0.8/ af, so the particle config- 
urations are severely jammed. We will measure space 
and time in units of o\ and To = (miaf/e) 1 ^ 2 . The tem- 
perature T will be measured in units of e/fcs, and the 
viscosity rj in units of ero / c 3 . Very long annealing times 
(~ 2.5 x 10 5 ) are chosen in our case. Then, for T > 0.267 
no appreciable aging effect is detected in various quanti- 
ties such as the pressure or the density time correlation 
function, whereas at the lowest temperature, T = 0.234, 
a small aging effect remains in the density time correla- 
tion function. 

Let us consider the incoherent density correlation 
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FIG. 1. rj/T versus r a at various temperatures. Those 
are obtained from nonequilibrium MD in shear flow |J. The 
straight line represents rj/T = 2nT a . 
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FIG. 2. Dr a versus r a . The solid line represents the 
Stokes-Einstein value DesTc = (2tv)~ 2 arising from Eq.(l). 



function, F s (q,t) = (X)j=i ex P [*Q ' Ar^ (t)]) /A^i for the 
particle species 1, where Arj(t) = Tj(t) — t\,(0) is the 
displacement vector of the j-th particle. This function 
may be introduced also in shear flow H. The a re- 
laxation time r a is then defined by F s (q,r a ) = e _1 
at q = 2ir for various T (and the shear rate 7). We 
also calculate the coherent time correlation function, 
Sn(q,t) = (m(q,t)ni(—q,0)), for the density fluctu- 
ations of the particle species 1. Interestingly, the de- 
cay profiles of Sn(q,t) at its first peak wave number 
q = Qm ~ 2-7T and F s (q, t) at q = 2ir nearly coincide 
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FIG. 3. 47rr- G s (r,t) versus r &t t — r a . The solid line is 
for T = 0.267 and the broken line is for T = 0.473. The 
dotted line represents the Brownian motion result. Note that 
the areas below the curves give 6Dr a . 

in the whole time region studied (t <2x 10 5 ) within 5%. 
Hence Sn(q mi T a )/ Sn(q m ,0) = e _1 holds for any T in 
our simulation. Such agreement is not obtained for other 
wave numbers, however. Furthermore, some neutron- 
spin-echo experiments JlTj showed that the decay time 
of Sn(q m ,t) is nearly equal to the stress relaxation time 
and as a result the viscosity 77 is of order r a . In accord 
with this experimental result, we obtain a simple linear 
relation in our simulation, 



(2na 1 /ql l )r]/kBT 



(1) 



in the original units. Fig.l shows that Eq.(l) is valid for 
any T and 7 over a wide range of r a . Here we may de- 
fine a g-dependent relaxation time T q by F(q,T g ) = e _1 . 
Thus, at the peak wave number q = q m , the effective 
diffusion constant D q = l/q 2 T q is given by the Einstein- 
Stokes form even in highly supercooled liquids. 

However, notice that the usual diffusion constant is 
the long wavelength limit, D = lim q ^Q D q . It is 
usually calculated from the mean square displacement, 
((Ar(i)) 2 ) = (ESi^iCO) 2 )/^!- The crossover of this 
quantity from the plateau behavior arising from motions 
in transient cages to the diffusion behavior 6Dt has been 
found to take place around t ~ 0.1t q Q. In Fig. 2 we 
plot Dr a versus r a , which clearly indicates breakdown of 
the Einstein-Stokes relation in agreement with the exper- 
imental trend. To examine the diffusion process in more 
detail we here introduce the van Hove self-correlation 
function, G s {r,t) = (£f^ S(A rj (t) - r))/JVi. Then, 



F*(q,t) 



qr 



(2) 



is the 3D Fourier transformation of G s (r,t). At q = 2tt, 
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FIG. 4. Mobile particles of the species 1 at T — 0.267. The time interval t is 0.125r Q in (a), r a in (b), and 10r a in (c). 
The radii of the spheres are |Arj(t)|A/ ((Ar(i)) 2 ) and the centers are at \[rj(ta) + rj(to + t)]. The system linear dimension 
is L — 23.2. The darkness of the spheres represents the depth in the 3D space. The heterogeneity is significant at 0.125tq, and 
T a but is much decreased at 10r Q . 



the integrand in Eq.(2) vanishes at r = 1 and the integral 
in the region r < 1 is confirmed to dominantly determine 
the decay of F s (27r, t) or r a . On the other hand, the mean 
square displacement 

poo 

((Ar(t)) 2 ) = / dr4nr 4 G s (r,t) (3) 
Jo 

is determined by the particle motions out of the cages. 
In Fig. 3 we display 4tit G s (r, r a ) versus r at zero shear, 
where t q = 3.2 and 2000 for T = 0.473 and 0.267, respec- 
tively. These curves may be compared with the Gaussian 
(Brownian motion) result, {2/ir) 1 / 2 £~ 3 r i exp(— r 2 /2£ 2 ), 
where 3£ 2 = QDEs T a = 3/2-7T 2 is the Einstein- Stokes 
mean square displacement. Because the areas below the 
curves give 6-Dt q , we recognize that the particle motions 
over large distances r > 1 are much enhanced at low T, 
leading to the violation of the Einstein-Stokes relation. 

We then visualize the heterogeneity of the difFusivity. 
To this end, we pick up mobile particles of the species 
1 with |Arj(£))| > £ c in a time interval [to, to + t] and 
number them as j = 1, • ■ • , N rn . Here £ c is defined such 
that the sum of Arj(t) 2 of the mobile particles is 66% of 
the total sum (= QDtN x for t > 0.1r Q ). In Fig.4 these 
particles are written as spheres with radius 

aj (t)^\ArM/Vi(^W) (4) 

located at Rj(t) = \[rj{to) + rj(t + t)] in three time 
intervals, [to, to + 0.125r Q ] in (a), [to, to + i~ a ] in (b), 
and [to, to + 10r Q ] in (c). The lower cut-off £ c and the 
mobile particle number N m increase with increasing t 
as 0.48 and 571 in (a), 1.3 and 725 in (b), and 2.9 
and 1316 in (c), respectively. Here they approach the 
Gaussian results, £ c = 0.403(t/r Q ) 1 / 2 and N m = 1800, 
for t 3> r a . In our case the ratio of the second mo- 
ments C2 = J2f=i a j(t) 2 / a j(t) 2 i s ne ld fixed at 0.66 
independently of t, while the ratio of fourth moments 
C4 = Ylf=i a j(t) 4 /J2f=i a j{i) i turns out to be close to 1 



as c 4 = 0.97 in (a), 0.92 in (b), and 0.90 in (c). We can 
see that the large scale heterogeneities in (a) and (b) are 
much weakened in (c) and that the areas of the spheres 
have the largest variance in (a). In fact, the variance de- 

fined by V = N m £^ % (t) 4 /(Ef=i ^M 2 ) 2 - 1 is 0.94 
in (a), 0.41 in (b), and 0.32 in (c). Here c 4 ->■ 0.833 
and V — > 0.13 for t » r a . The above result is con- 
sistent with the fact that the non-Gaussian parameter 
A2(t) takes a maximum of 3.1 at t = 0.125t q at this 
temperature. This is because the statistical average of 
V (taken over many initial times to) is related to ^2(t) 
by (V) £* (5(c 4 )(iV m )/3c 2 .7V 1 )(l + A 2 (t)) - 1. We may 
also conclude that the significant rise of A 2 {t) in glassy 
states originates from the heterogeneity in accord with 
some experimental interpretations [Q. 

Furthermore, we consider the Fourier component of the 
diffusivity density defined by 

V q (to,t) = J2 aj (t) 2 exp[iq ■ (r - Rj{t))}, (5) 

7=1 

which depends on the initial time to and the final 
time to + t. The correlation function Sx>{q,t,T) = 
(f q(to + T, t)T>-q(to, t)) is then obtained after averag- 
ing over many initial states. We confirm that S-r>(q,t,T) 
tends to its long wavelength limit for q < where £ 
coincides with the correlation length of the heterogene- 
ity structure of the bond breakage ||[l2]]. As the differ- 
ence t of the initial times increases with fixed t = r Q , 
Sx>{q,T a ,T) relaxes as exp[— (r/r^) c ] for q < where 
c ~ 0.5 at T = 0.267 and t/j ~ 3r a is the life time of the 
heterogeneity structure. The two-time correlation func- 
tion among the broken bond density P^,pl also relaxes 
with Th in the same manner. 

We naturally expect that the distribution of the par- 
ticle displacement Arj (t) in the active regions should be 
characterized by the local diffusion constant D(x,t) de- 
pendent on the spatial position x = (x, y, z) and the 
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FIG. 5. A test of the scaling plot ^6Dt4-Kr 2 G s (r,t) ver- 
sus r/V&Dt for 0.5r a < t < 3r a (t = (0.5 + 0.25n)r Q , 
n = 0, 1, 2, • • • , 10). The inset shows the curves at longer 
times, 10r a < t < 100r Q (t = 10nr a , n = 1,2, • • • , 10), where 
the heterogeneity effect is smoothed out. The dotted lines are 
the Gaussian form. 

time interval t. The van Hove correlation function 
G s (r, t) may then be expressed as the spatial aver- 
age of a local function G s (x,r,t), which is given by 
[AnD(x,t)t]- 3 / 2 exp[-r 2 /AD(x,t)t] in the relatively ac- 
tive regions. To check this conjecture we plot the scaled 
function V 6Dt4nr 2 G s (r, t) versus r* — r/>/6Dt in Fig. 5. 
The areas below the curves are fixed at 1. At relatively 
short times t < 3r a , the curves in the region r > 1 or 
r* > {QDt)- 1 / 2 , which give dominant contributions to 
((Ar(t)) 2 ), tend to a master curve quite different from 
the rapidly decaying Gaussian tail. Note that the peak 
position of each curve at larger r* corresponds to r = 1 in 
Fig. 5. This asymptotic law is consistent with the picture 
of the space-dependent diffusion constant in the active 
regions. It is also important that the heterogeneity struc- 
ture remains unchanged in the time region t < 77, ~ 3r Q . 
At longer times t > 10r Q the curves approach the Gaus- 
sian form as can be seen in the inset of Fig. 5. Of course, 
4irr 2 G s (r, i) for r < 1 does not scale in the above manner, 
because it is the probability density of a tagged particle 
staying within a cage. This short-range behavior deter- 
mines the decay of F s (2ir,t) as noted below Eq.(2). 
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